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Program Summary: 

• Title of program: KCpKr], Version 1.0 (July 1996) 

• Program obtainable by anonymous ftp from the host 

crunch. ikp.physik.th-darmstadt.de in the directory pub/ohl/circe. 

• Licensing provisions: Free software under the GNU General Public 
License. 

• Programming language used: Fortran?? 

• Number of program lines in distributed prograun, including test 
data, etc.: w 1100 (excluding comments) 

• Computer/Operating System: Any with a Fortran?? programming 
environment. 

• Memory required to execute with typical data: Negligible on the 
scale of typical applications calling the library. 

• Typical running time: A small fraction (typically a few percent) of the 
running time of applications calling the library. 

• Purpose of program: Provide simple and reproducible, yet realistic, 

parameterizations of the e*- and 7-beam spectra for linear colliders. 

• Nature of physical problem: The intricate beam dynamics in the inter- 
action region of a high luminosity linear collider at ^/s = 500GeV result in 
non-trivial energy spectra of the scattering electrons, positrons and pho- 
tons. Physics simulations require simple and reproducible, yet realistic, 

parameterizations of these spectra. 

• Method of solution: Parameterization, curve fitting, Monte Carlo event 
generation. 

• Keywords: Event generation, beamstrahlung, linear colliders. 
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1 Introduction 



Despite the enormous quantitative success of the electro-weak standard model 
up to energies of 200GeV, neither the nature of electro-weak symmetry break- 
ing (EWSB) nor the origin of mass are understood. 

From theoretical considerations, we know that clues to the answer of these 
open questions are hidden in the energy range below Aewsb — 47rw « S.lTeV. 
Either we will discover a Higgs particle in this energy range or signatures for 
a strongly interacting EWSB sector will be found. Experiments at CERN's 
Large Hadron Collider (LHC) will shed a first light on this regime in the next 
decade. In the past is has been very fruitful to complement experiments at 
high energy hadron colliders with experiments at e+e^-colliders. The simpler 
initial state allows more precise measurements with smaller theoretical errors. 
Lucid expositions of the physics opportunities of high energy e"*"e~ colliders with 
references to the literature can be found in 

However, the power emitted by circular storage rings in form of synchrotron 
radiation scales like {E/mY/R^ with the energy and mass of the particle and 
the radius of the ring. This cost becomes prohibitive after LEP2 and a Linear 
Collider (LC) has to be built instead. 

Unfortunately, the "interesting" hard cross sections scale like 1/s with the 
square of the center of mass energy and a LC will have to operate at extremely 
high luminosities in excess of 10^'^cm^^s^^. To achieve such luminosities, the 
bunches of electrons and positrons have to be very dense. Under these con- 
ditions, the electrons undergo acceleration from strong electromagnetic forces 
from the positron bunch (and vice versa). The resulting synchrotron radiation is 
called beamstrahlung and has a strong effect on the energy spectrum D{xi, X2) 
of the colliding particles. This changes the observable e+e~ cross sections 

° "(*) ^ — IFi—i'^)^ dxidx2De+^-{xi,X2;Vs)J{n\V,)^——{xiX2s) 



(la) 



and produces luminosity for e^7 and 77 collisions: 

(g) = / dxidx2D^±^{xi,X2]^)J{n' ,n)—^^{xiX2s) 

(lb) 

——(^s)^ I dxidx2Dj^{xi,X2;Vs)J{n\n)—^{xiX2s) 

(Ic) 

Therefore, simulations of the physics expected at a LC need to know the spectra 
of the and 7 beams precisely. 

Microscopic simulations of the beam dynamics are available (e.g. ABEl|^, 
CAIN[Q and Guinea-Pig|^) and their predictions are compatible with each 
other. But they require too much computer time and memory for direct use 
in physics programs. Kt'pKry provides a fast and simple parameterization of the 
results from these simulations. Furthermore, even if the computational cost of 
the simulations would be negligible, the input parameters for microscopic sim- 
ulations are not convenient for particle physics applications. Due to the highly 
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Table 1: Accelerator parameters for three typical designs at ^fs — 500GeV 
and \fs = ITeV. The resulting distributions are shown in figure ^. The design 
efforts are currently concentrated on a 350GeV-800GeV LC. Therefore the Tesla 
parameters for ITeV are slightly out of date. 



non-linear beam dynamics, the optimization of LC designs is a subtle art 
that is best practiced by the experts. Furthermore, particle physics applications 
need benchmarking and easily reproducible parameterizations are required for 
this purpose. 

The parameterizations in KcpKi] are not based on approximate solutions (cf . ) 
of the beamstrahlung dynamics. Instead, they provide a "phenomenological" 
description of the results from full simulations. The parameterizations are as 
simple as possible while remaining consistent with basic physical principles: 

1. positivity: the distribution functions D{xi,X2) must not be negative in 
the physical region [0, 1] x [0, 1]. 

2. integrability: the definite integral of the distribution functions over the 
physical region [0, 1] x [0, 1] must exist, even though the distributions can 
have singularities. 

This paper is organized as follows: I start in section 2 with a discussion of 
the input for the microscopic simulations. In section 3 I describe the usage 
of the Klpkt] library and in section || I discuss some technical details of the 
implementation. After discussing the parameterizations available in version 1.0 
in section ||, 1 conclude in section |6[ 

2 Parameters 

The microscopic simulation program Guinea-Pig used for the current version 
of the parameterizations in Y^ipmrj simulates the passage of electrons through a 



4 



= 500GeV ^ = SOOGeV 




1 - Xe± 



Figure 1: Version 1, revision 1996 07 11 of the factorized e^- and 7-distributions 
at ^/s = 500GeV and ^/s = ITeV in a doubly logarithmic plot. The accelerator 
parameters are taken from table nl 
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Table 2: Accelerator parameters for the Tesla design at three planned [g ener- 
gies. The resulting distributions are shown in figure 




Figure 2: Version 1, revision 1996 07 11 of the factorized e - and 7-distributions 
for Tesla in a doubly logarithmic plot. The accelerator parameters are taken 
from table ^ 
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bunch of electrons (and vice versa) . It takes the foUowing accelerator parameters 
as input: 

E : the energy of the particles before the beam-beam interaction. 

-^particles '■ the number of particles per bunch. 

ex,y '■ the normalized horizontal and vertical emittances. 

(3* y : the horizontal and vertical beta functions. 

'^x.y,z '■ the horizontal, vertical and longitudinal beam size. A Gaussian shape 
is used for the charge distribution in the bunches. 

/rep ■ the repetition rate. 

"•bunch ■ the number of bunches per train. 

The transversal beam sizes, beta functions and normalized emittances for rela- 
tivistic particles are related by 

Ply-^— (2) 

The parameters used in the most recent revision of the parameterizations are 
collected in tables |] and ||. The resulting factorized electron/positron and pho- 
ton distributions in version 1 of the parameterizations are depicted in figures ^ 
and 1^. 

The most important purpose of KipKi] is to map the manifold of possible 
beam spectra for the NLC to a finite number of reproducible parameterizations. 
The distributions 

D'^Z(''i,X2;V~s) (3) 
provided by KcpKri are indexed by three integers 

a : the accelerator design class: currently there are three options: S-band 

Tesla 1^, X-band |ll|. More variety will be added later, in particular 
the e~e^ mode and the e^7 and 77 laser backscattering modes of these 
designs. 

ly : the version of the parameterization: over the years, the form of the param- 
eterizations can change, either because better approximations are found 
or because new simulation programs become available. All versions will 
remain available in order to be able to reproduce calculations. 

p : the revision date for the parameterization: a particular parameterization 
can contain bugs, which will be fixed in subsequent revisions. While only 
the most recent revision should be used for new calculations, old revisions 
will remain available in order to be able to reproduce calculations. 

The continuous parameter -y/s in (^) is misleading, because accelerator param- 
eters have been optimized for discrete values of the energy. Therefore the dis- 
tributions are not available for all values of ^/s. 

The usage of the distributions in application programs is discussed in sec- 



tion 3.1. KbpnT] provides for each of the distributions a non- uniform random 



variate generator, that generates energy fractions according to the distributions. 



The usage of these generators is discussed in section 3.2 
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3 Usage 



3.1 Distributions 

A generic interface to all distributions Dp-^p^(xi, X2) is given by the circe func- 
tion 

{API documentation) = 

double precision circe, d, xl, x2 

integer pi, p2 

d = circe (xl , x2, pi, p2) 

where the energy fractions are specified by xi^2 and the particles pi,2 are iden- 
tified by their standard Monte Carlo codes :[p^ 

{Particle codes) = 

integer ELECTR, POSITR, PHOTON 
parameter (ELECTR = 11) 
parameter (POSITR = -11) 
parameter (PHOTON = 22) 

The distributions can have integrable singularities at the end points, therefore 
the calling functions must not evaluate them at the endpoints and 1. This 
is usually not a problem, since standard mapping techniques (cf. ( |l0[) below) 
will have to be used to take care of the singularity anyway. Nevertheless, all 
applications should favor open quadrature formulae (i.e. formulae not involving 
the endpoints) over closed formulae. The distributions are guaranteed to vanish 
unless < xi^2 < 1, with two exceptions. Firstly, the value —1 allows to pick 
up the integral of the continuum contribution: 

Dpip2{-^,X2) = lim / dxi Dp^p2{xi,X2) (4a) 

Dpip2{xi,-l) ^ / dx2Dp^p2{xi,X2) (4b) 

-OpiP2(-l) -1) = 1™ / dxidx2 Dp^p2{xi,X2) (4c) 

The other exception is that the strength of (5-function contributions at the end- 
point can be picked up from the value at this endpoint: 

Z^e+e- {xi,X2) — -De+e- (1, 1)S{1 — Xi)5{l — ^2) + smooth and single S (5a) 
De±-f{xi, X2) = Df.±^{l, X2)5{\ — Xl) + smooth (5b) 
D^^± {xi, X2) = D-ye^ i^ii 1)<5(1 ^ ^2) + smooth (5c) 

The use of these special values is demonstrated in an example in section |3.1.l| 
below. 

The distributions are normalized such that 

i-l+e 

lim / dxidx2 Dg+^- {xi, X2) = 1- (6) 
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and the nominal e"'"e~-luminosity of the currently active accelerator design can 
be retrieved from the database with the subroutine circel. The value is given 
in units of 

Vo-^y-^ = lO^^^cm-^sec"^ (7) 

where v = lO^'sec « year/7r is an "effective year" of running with about 30% 
up-time . 

{API documentation) += 

double precision lumi 
call circel (lumi) 

A particular parameterization is selected by the circes function: 

{API documentation) += 

double precision xlm, x2m, roots 
integer acc, ver, rev, chat 

call circes (xlm, x2m, roots, acc, ver, rev, chat) 

The parameter roots corresponds to the nominal center of mass energy ^/s / GeV 
of the colhder. Currently = 350GeV, SOOGeV, SOOGeV, ITeV (i.e. 350D0, 
500D0, 800D0 and lOOODO) are supported. Application programs can not assume 
that energy values are interpolated. For convenience, e.g. in top threshold scans 
around 350GeV, a small interval around the supported values will be accepted 
as synonymous with the central value, but a warning will be printed. Section ^ 
should be consulted for the discrete values supported by a particular version of 
the parameterizations. Negative values of roots will keep the currently active 
value for ^/s. 

The parameters xlm and x2m will set thresholds ^l.min cLIld 2^2, min 

for the 



event generation in the routines described in section 3.2. 

The parameter acc selects the accelerator design. Currently the following 
accelerator codes are recognized: 

{Accelerator codes) = 

integer SBAND, TESLA, XBAND 

parameter (SBAND = 1, TESLA = 2, XBAND = 3) 
integer NACC 
parameter (NACC = 3) 

Negative values will keep the currently active accelerator. Later I will add 
the e~e~ mode and the e^7 and 77 laser backscattering modes of these designs: 

{Future API documentation) = 

integer SBAND, TESLA, XBAND 

parameter (SBAND = 1, TESLA = 2, XBAND = 3) 

integer SBNDEE, TESLEE, XBNDEE 

parameter (SBNDEE = 4, TESLEE = 5, XBNDEE = 6) 

integer SBNDEG, TESLEG, XBNDEG 

parameter (SBNDEG = 7, TESLEG = 8, XBNDEG = 9) 

integer SBNDGG, TESLGG, XBNDGG 

parameter (SBNDGG = 10, TESLGG =11, XBNDGG = 12) 
integer NACC 
parameter (NACC = 12) 
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The ver parameter is used to determine the version as follows: 



ver > : a frozen version which is documented in section ^. For example, 
version 1 is a family of factorized Beta distributions: D{xi, X2) oc a;°^ (1 — 

Xi)^^X^Hl-X2)'''. 

ver — : the latest experimental version, which is usually not documented and 
can change at any time without announcement. 

ver < : keep the currently active version. 

The rev parameter is used to determine the revision of a version as follows: 

rev > : a frozen revision which is documented in section |^. The integer rev is 
constructed from the date as follows: rev = 10'* • year + 10^ • month + day, 
where the year is greater than 1995. Since Fortran?? ignores whitespace, 
it can be written like 1996 07 11 for readability. If there is no exact 
match, the most recent revision before the specified date is chosen. 

rev = : the most recent revision. 

rev < ; keep the currently active revision. 

Finally, the parameter chat controls the "chattiness" of circe. If it is 0, only 
error messages are printed. If it is 1, the parameters in use are printed whenever 
they change. Higher values of chat can produce even more diagnostics. 

In addition to the generic interface circe, there are specialized functions for 
particular particle distributions. Obviously 

-De±7 {X1,X2,S) ^ D°!^ll (X2 , , s) (8) 

and there are three independent functions I?e-e+, D^-^ and D^^ for the e^e~ 
colliders with reasonable mnemonics: 

{API documentation) += 

double precision circee, circeg, circgg 
d = circee (xl, x2) 
d = circeg (xl, x2) 
d = circgg (xl, x2) 

Calling the latter three functions is marginally faster in the current implemen- 
tation, but this can change in the future. 



3.1.1 Example 

For clarification, let me give a simple example. Imagine we want to calculate 
the integrated production cross section 

<^x{s) ^ J dXidX2Cre+e--,x{xiX2s)D^+^-{xi,X2,s) (9) 

Since the distributions are singular in the Xi^2 ~* 1 limit, we have to map away 
this singularity with 

x-^^ = (l-a;)!/" (10a) 



10 



Therefore 




(10b) 



with T] sufficiently large to give the integrand a finite limit at a; — > 1. If / 
diverges like a power f{x) oc 1/(1 — x)^ , this means 77 > 1/(1 — /3). 

As a specific example, let us "measure" a one particle s-channel exchange 
cross section 



double precision function sigma (s) 

implicit none 

double precision s 

sigma = IdO / s 

end 



I will present the example code in a bottom-up fashion, which should be intuitive 
and is described in some more detail in appendix Assuming the existence of a 
one- and a two-dimensional Gaussian integration function gaussl and gauss2j^ 
we can perform the integral as follows: 

{Gauss integration) = 

s = sigma (IdO) * circee (IdO, IdO) 
$ + gaussl (dl, OdO, IdO, EPS) 
$ + gaussl (d2, OdO, IdO, EPS) 
$ + gauss2 (dl2, OdO, IdO, OdO, IdO, EPS) 
write (*, 1000) ' delta(sigma) (Gauss) =', (s-ldO) *100dO 
1000 format (IX, A22, IX, F5.2, '7.') 

Note how the four combinations of continuum and (5-peak are integrated sep- 
arately, where you have to use three auxiliary functions dl, d2 and dl2. The 
continuum contribution, including the Jacobian: 

(sEunple . f )-|-= 

double precision function dl2 (tl, t2) 
implicit none 

double precision tl, t2, xl, x2, sigma, circee 
(EPS & PWR ) 
xl = IdO - tl**PWR 
x2 = IdO - t2**PWR 
dl2 = PWR*PWR * (tl*t2)**(PWR-ldO) 
$ * sigma (xl*x2) * circee (xl, x2) 

end 

the first product of continuum and (5-peak: 
(scunple . f )+= 



cr(s) (X - 



(11) 



(sample . f ) 



double precision function dl (tl) 
implicit none 

double precision tl, xl , sigma, circee 



^They are provided in the example program sample. f. 
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(EPS & PWR ) 

xl = IdO - tl**PWR 

dl = PWR * tl**(PWR-ldO) * Sigma (xl) * circee (xl, IdO) 
end 

and the second one: 

(sELmple . f )+= 

double precision function d2 (t2) 
implicit none 

double precision t2, x2, sigma, circee 

(EPS & PWR ) 

x2 = IdO - t2**PWR 

d2 = PWR * t2**(PWR-ldO) * sigma (x2) * circee (IdO, x2) 
end 

Below you will see that the power of the singularity of the e+e^ distributions 
at x ^ 1 is w —2/3. To be on the safe side, we choose the power 77 in ( p^ 
as 5. It is kept in the parameter PWR, while EPS is the desired accuracy of the 
Gaussian integration: 

(EPS & PWR) = 

double precision EPS, PWR 
parameter (EPS = ld-6, PWR = 5d0) 

These code fragments can now be used in a main program that loops over 
energies and accelerator designs 

(sample . f )+= 

program sample 
implicit none 

{Accelerator codes ) 
(EPS & PWR ) 

{Other variables in seunple ) 
integer acc, ver, i 
double precision roots (2) 
data roots / 500D0, lOOODO / 
do 10 acc = 1, NACC 
do 11 ver = 1, 1 
do 12 i = 1, 2 

call Circes (OdO, OdO, roots(i), acc, ver, 1996 07 29, 
{Gauss integration ) 
{Monte Carlo integration ) 
12 continue 
11 continue 
10 continue 
end 

with the following result 

{Sample output)= 

circe :message : starting up . . . 

Circe :message: Id: circe. nw,v 1.22 1996/07/27 19:52:28 ohl Exp 
circe :message : updating 'acc' to SBAND 
circe :message : updating 'ver' to 1 
circe :message : updating 'rev' to 19960729 



12 



deltaCsigma) (Gauss) = 3.797o 

delta(sigma) (MC) = 3.747. 

+/- . 067. 

circe :message : updating 'roots' to 1000.0 
delta(sigma) (Gauss) = 10.117. 

delta (sigma) (MC) = 9.977. 

+/- . 147. 

circe : message : updating 'roots' to 500.0 

circe : message : updating 'acc' to TESLA 

delta(sigma) (Gauss) = 3.117. 

delta (sigma) (MC) = 3.097. 

+/- .047. 

circe :message : updating 'roots' to 1000.0 

delta(sigma) (Gauss) = 3.987. 

delta(sigma) (MC) = 3.967. 

+/- .077. 

circe :message : updating 'roots' to 500.0 

circe : message : updating 'acc' to XBAND 

delta(sigma) (Gauss) = 4.967. 

delta (sigma) (MC) = 4.967. 

+/- . 107. 

circe :message : updating 'roots' to 1000.0 
delta(sigma) (Gauss) = 21.317. 
delta(sigma) (MC) = 21.727. 

+/- .457. 

We almost forgot to declare the variables in the main program 

{Other variables in sample)= 
double precision s 

double precision gaussl, gauss2, circee, sigma, dl, d2, dl2 
external dl, d2, dl2 

This concludes the integration example. It should have made it obvious how to 
proceed in a realistic application. 

In section 3.2.1 below, I will describe a Monte Carlo method for calculating 
such integrals efficiently. 



3.2 Generators 

The function circe and its companions are opaque to the user. Since they 
will in general contain singularities, applications will not be able to generate 
corresponding samples of random numbers efhciently To fill this gap, four 
random number generators are provided. The subroutine girce will generate 
particle types pi,2 and energy fractions xi^2 in one step, according to the selected 
distribution. Particle pi will be either a positron or a photon and p2 will be 
either an electron or a photon. The energy fractions are guaranteed to be above 
the currently active thresholds: Xi > Xi^^nin- This can be used to cut on soft 
events — the photon distributions are rather soft — which might not be interesting 

^The implementation of the flavor selection with non- vanishing thresholds xi 
and X2 min is moderately inefliicient at the moment. It can be improved by a factor of two. 
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in most simulations. 

{API documentation) += 

call girce (xl, x2, pi, p2, rng) 

The output parameters of girce are identical to the input parameters of circe, 
with the exception of rng. The latter is a subroutine with a single double 
precision argument, which will be assigned a uniform deviate from the interval 
[0, 1] after each call: 

{API documentation) += 

subroutine rng (r) 
double precision r 

r = {uniform deviate on [0, 1] ) 
end 

Typically, it will be just a wrapper around the standard random number gener- 
ator of the application program. For studies with a definite initial state, three 
generator functions are available. 

{API documentation) 

call gircee (xl, x2, rng) 

call girceg (xl, x2, rng) 

call gircgg (xl, x2, rng) 



3.2.1 Example 



Returning to the example from section 3.2.1 , I present a concise Monte Carlo 
algorithm for calculating the same integral: 

{Monte Carlo integration)— 
s = OdO 
s2 = OdO 

do 100 n = 1, NEVENT 

call gircee (xl, x2, random) 

w = sigma (xl*x2) 

s = s + w 

s2 = s2 + w*w 
100 continue 

s = s / dble (NEVENT) 
s2 = s2 / dble (NEVENT) 

write (*, 1000) 'delta(sigma) (MC) =', (s-ldO) *100dO 
write (*, 1000) ' +/-', 
$ sqrt ( (s2-s*s) /dble (NEVENT) ) *100d0 

{Other variables in sample)+= 

double precision w, s2, xl, x2 
external random 
integer NEVENT, n 
parameter (NEVENT = 10000) 

Here is a simple linear congruential random number generator for the sample 
program. Real applications will use their more sophisticated generators instead. 

(sample . f )+= 

subroutine random (r) 
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implicit none 
double precision r 
integer m, a, c 

parameter (M = 259200, A = 7141, C = 54773) 

integer n 

save n 

data n /O/ 

n = mod(n*a+c,m) 

r = dble (n) / dble (m) 

end 

If the cross section is slowly varying on the range where the xi^2 distributions 
are non-zero, this algorithm is very efficient. 

However, if this condition is not met, the explicit form of the parametcri- 
zations in section ^ should be consulted and appropriate mapping techniques 
should be applied. The typical example for this problem is a narrow resonance 
just below the nominal beam energy. 



3.2.2 Event Generators 

For Monte Carlo event generators that use the standard /hepevt/ common 
block [Q, the addition of the KipKr] library is trivial. During the initialization 
of the event generator, the circes subroutine is called to set up KtpKjy's internal 
state. For example: 

(Initialize event generator) = 

call circes (OdO, OdO, roots, acc, ver, 1996 07 11, 1) 

During event generation, before setting up the e^^e^ initial state, the gircee 
subroutine is called with the event generator's random number generator: 

{Event generation) = 

call gircee (xl, x2, random) 

The resulting energy fractions xi and X2 are now available for defining the initial 
state electron 

{Event generation) += 

isthep(l) = 101 
idhep(l) = ELECTR 
p}iep(l,l) = OdO 
phep(2,l) = OdO 
phep(3,l) = xl * ebecun 
pliep(4,l) = xl * ebeami 
phep(5,l) = OdO 

and positron. 

{Event generation) += 

istliep(2) = 102 
idliep(2) = POSITR 

pliep(l,2) = OdO 

pliep(2,2) = OdO 

phep(3,2) = - x2 * ebeam 

phep(4,2) = x2 * ebesun 

pliep(5,2) = OdO 
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data statements in 
circes, forming a 
database of 
parameterizations. 



distribution 

D{xi,X2) 
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/circec/ common 
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preloaded 
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random 
variables: 

Xi,X2 



Figure 3: Architecture of Klpkt]: circes () selects energy and accelerator and 
loads the parameterization. The function circeO calculates the values of the 
selected distribution function at the given energy fractions. The subroutine 
girceO generates energy fractions using a specified random number generator 
in accordance with the selected distribution. 



Using KipKij with other event generators should be straightforward as well. 



4 Technical Notes 

The structure of Kipur] is extremely simple (cf. figure ^ and is mainly a book- 
keeping excercise. All that needs to be done is to maintain a database of avail- 
able parameterizations and to evaluate the corresponding functions. The only 
non trivial algorithms are used for the efficient generation of random deviates. 

I have avoided the use of initialized common blocks (i.e. block data sub- 
routines), because the Fortran?? standard does not provide a portable way of 
ensuring that block data subroutines are actually executed at loading time. 
Instead, the /circom/ common block is tagged by a "magic number" to check 
for initialization and its members are filled by the circes subroutine when 
necessary. 

A more flexible method would be to replace the data statements by reading 
external files. This option causes portability problems, however, because I would 
have to make sure that the names of the external files are valid in all files systems 
of the target operating systems. More significantly, splitting the implementation 
into several parts forces the user to keep all files up to date. This can be a 
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Table 3: Version 1, revision 1996 07 11 of the beam spectra at 500 GeV. The rows 
correspond to the luminosity per effective year, the integral over the continuum 
and the powers in the factorized Beta distributions (E2|). 
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Table 4: Version 1, revision 1996 07 29 of the beam spectra at 1 TeV. 

problem, because Fortran source files and data input files will typically be kept 
in different parts of the file system. 

The option of implementing Kt'pKjy statelessly, i.e. with pure function calls 
and without common blocks, has been dismissed. While it would have been 
more straightforward on the side of the library, it would have placed the burdon 
of maintaining state (accelerator, energy, etc.) on the application program, 
thereby complicating them considerably. Keeping an explicit state in KCpKi] 
has the additional benefit of allowing to precompute certain internal variables, 
resulting in a more efficient implementation. 

5 Parameterizations 

Version 1.0 of KtpKry supports just one version of the parameterizations. Future 
versions will provide additional parameterizations. 
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Table 5: Version 1, revision 1996 07 29 of the beam spectra for TESLA. 




1 0.01 10-4 10-6 10-8 10-8 10-6 10-4 0.01 1 



1 - Xe± x^ 

Figure 4: Fit of the e^- and 7-distributions for the S-Band design at y/s = 
500GcV. The open circles with error bars are the result of the Guinea-Pig 
similation. The full line is the fit. 
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Figure 5: Fit of the e^- and 7-distributions for the Tesla design at ^/s = 
500GeV. 




Figure 6: Fit of the e^- and 7-distributions for the X-Band design at ^/s = 
500GeV. 
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Tesla, = ITeV 



Tesla, y/s = ITeV 




Figure 7: Fit of the e^- and 7-distributions for the Tesla design at -^i = ITcV. 
5.1 Version 1 

The first version of the parameterization uses a simple factorized ansatz 

D'^lP^{xuX2. s) = df/{x,)df/{x2) (12a) 

where the distributions are simple Beta distributions: 

(P^lPix) = ao''(5(l -x)+ a^x^" (1 - x)<' (12b) 
d^ip(x) = a^x^" (1 - xY"^' (12c) 

This form of the distributions is motivated by the observation Q that the dis- 
tributions diverge like a power for x 1 and vanish at a: — > 0. The behavior of 
the 7 distributions is similar with the borders exchanged. 

5.1.1 Fitting 

The parameters in ( |l^ have been obtained by a least-square fit of (|l|) to 
histograms of simulation results from Guinea-Pig. Some care has to taken when 
fitting singular distributions to histogrammed data. Obviously equidistant bins 
are not a good idea, because most bins will be almost empty (cf. figures 
and H) and consequently a lot of information will be wasted. One solution to 
this problem is the use of logarithmic bins. This, however, maps the compact 
region [0,1] x [0,1] to [— cx),0] x [— oo,0], which is inconvenient because of the 
missing lower bounds. 

The more appropriate solution is to use two maps 

: [0, 1] - [0, 1] 

X y — X ' ' 

where x — x.^ ov x — \ — Xf,± , and to bin the result equidistantly. If 77 is chosen 
properly (cf. (p^), the bin contents will then fall ofi^ at the singularity. The 
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fits in tables |j, and ^ have been performed with 77 = 5 and the resulting bin 
contents can be read off from figures |^-|7[ 

Using this procedure for binning the results of the simulations, the popu- 
lar fitting package MINUIT converges quickly in all cases considered. The 
resulting parameters are given in tables ^, ^, and ||. Plots of the correspond- 
ing distributions have been shown in figures and I It is obvious that an 
ansatz like ( [l^ is able to distinguish among the accelerator designs. Thus it 
can provide a solid basis for physics studies. 

In figures U-j?] I give a graphical impression of the quality of the fit, which ap- 
pears to be as good as one could reasonably expect for a simple ansatz like ( [l^ ) . 
Note that the histograms have non-equidistant bins and that the resulting Ja- 
cobians have not been removed. Therefore the bin contents falls off at the 
singularities, as discussed above. 

The errors used for the least-square fit had to be taken from a Monte 
Carlo (MC) study. Guinea-Pig only provides the -^/ri from Poissonian statistics 
for each bin, but the error accumulation during tracking the particles through 
phase space is not available. The MC studies shows that the latter error domi- 
nates the former, but appears to be reasonably Gaussian. A complete MC study 
of all parameter sets is computationally expensive (more than a week of pro- 
cessor time on a fast SGI). From an exemplary MC study of a few parameter 
sets, it appears that the errors can be described reasonably well by rescaling 
the Poissonian error in each bin with appropriate factors for electrons/positrons 
and photons and for continuum and delta. This procedure has been adopted. 

The x^/d.o.f.'s of the fits are less than 0(10). The simple ansatz (|l2|) is 
therefore very satisfactory. In fact, trying to improve the ad-hoc factorized 
Beta distributions by the better motivated approximations from ||^ or |l^, it 
turns out jl^ that ( p^ ) provides a significantly better fit of the results of the 
simulations. The price to pay is that the parameters in ( p^ ) have no direct 
physical interpretation. 



5.1.2 Generators 

For this version of the parameterizations we need a fast generator of Beta dis- 
tributions: 

(^''•''{x) o^x''-\l-x)''-^ (14) 

This problem has been studied extensively and we can use a published algo- 
rithm [|l^ that is guaranteed to be very fast for all a.b such that < a < 1 < 6, 
which turns out to be always the case (cf. tables I, I, and I). 



5.2 Future Versions 

There are two ways in which the parameterizations can be improved: 

more complicated functions: the factorized fits can only be improved marginally 
by adding more positive semi-definite factors to (p^). More improvement 
is possible by using sums of functions, but in this case, the best fits violate 
the positivity requirement and have to be discarded. 
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correlations: the parameterization in section 5.1 is factorized. While this 
is a good approximation, the simulations nevertheless show correlations 
among xi and X2- These correlations can be included in a future version. 



interpolation: the parameterization in section 5.1 is based on fitting the sim- 
ulation results by simple functions. Again, this appears to be a good 
approximation. But such fits can not uncover any fine structure of the 
distributions. Therefore it will be worthwhile to study interpolations of 
the simulation results in the future. A proper interpolation of results with 
statistical errors is however far from trivial: straightforward polynomial 
or spline interpolations will be oscillatory and violate the positivity re- 
quirement. Smoothing algorithms have to be investigated in depth before 
such a parameterization can be released. 

other simulations: besides js), other simulation codes are invited to con- 
tribute their results for inclusion in the KtpKjy library. 



6 Conclusions 

I have presented a library of simple parameterizations of realistic e^- and 7- 
beam spectra at future linear e''"e~-colliders. The library can be used for inte- 
gration and event generation. Emphasis is put on simplicity and reproducibility 
of the parameterizations for supporting reproducible physics simulations. 
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A Literate Programming 
A.l Paradigm 

I have presented the sample code in this paper using the literate programming 
paradigm. This paradigm has been introduced by Donald Knuth [|l9| and his 
programs TJ^X [|2^ and METRFONT [^ provide excellent examples of the virtues 
of literate programming. Knuth summarized his intention as follows ([^, p. 99) 

"Let us change our traditional attitude to the construction of pro- 
grams. Instead of imagining that our main task is to instruct a 
computer what to do, let us concentrate rather on explaining to 
human beings what we want a computer to do." 
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Usually, literate programming uses two utility programs to produce two kinds 
of files from the source 

tangle produces the computer program that is acceptable to an "illiterate" 
(Fortran, C, etc.) compiler. This process consists of stripping documen- 
tation and reordering code. Therefore it frees the author from having to 
present the code in the particular order enforced by a compiler for purely 
technical reasons. Instead, the author can present the code in the order 
that is most comprehensible. 

weave produces a documents that describes the program. Extensive cross refer- 
encing of the code sections is usually provided, which has been suppressed 
in this paper. If a powerful typesetting system (such a Tf^) is used, 
the document can present the algorithms in clear mathematical notation 
alongside the code. These features improve readability and maintainabil- 
ity of scientific code immensely. 



A. 2 Practice 

KtpKry uses the noweb system. This system has the advantage to work 
with any traditional programming language and support the essential features 



described in section A.l with minimal effort, noweb's tangle program only 
reorders the code sections, but does not reformat them. Therefore its output 
can be used just like any other "illiterate" program. 

The examples above should be almost self-explaining, but in order to avoid 
any ambiguities, I give another example: 
{Literate programming example) — 

{Code that has to be at the top ) 

{Other code ) 

I can start the presentation with the first line of the "other code" : 
{Other code) = 

line 1 of the other code 
If appropriate, the first line of the code that has to appear before the other code 
can be presented later: 
{Code that has to be at the top) = 

line 1 of the code at the top 
Now I can augment the sections: 
{Other code)+= 

line 2 of the other code 

{Code that has to be at the top)+= 

line 2 of the code at the top 
The complete "program" will be presented to the compiler as 

line 1 of the code at the top 
line 2 of the code at the top 
line 1 of the other code 
line 2 of the other code 

The examples in section |3.1.l| show that this reordering is particularly useful for 
declaring variables when they are first used (rather than at the beginning) and 
for zooming in on code inside of loops. 
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B Fortran Name Space 



In addition to the ten procedures and one common block discussed in section ^ 

• circe, circee, circeg, circgg, 

• girce, gircee, girceg, gircgg, 

• circes, circel, /circom/, 

there are two more globaUy visible functions which are used internally: 

• circem: error message handler, 

• girceb: efficient Beta distribution generator. 

Even if the /circom/ is globally visible, application programs must not ma- 
nipulate it directly. The circes, subroutine is provided for this purpose and 
updates some internal parameters as well. 

With features from the current Fortran standard (Fortran90), I could have 
kept the last two functions and the common block private. But since FortranOO 
has only been adopted by a small fraction of the high energy physics community, 
I have decided to remain in the confines of Fortran?? (except for the ubiquitous 
implicit none). 

Application programs wishing to remain compatible with future versions of 
KtpKjy must not use common blocks or procedures starting with circe or girce. 

C Updates 

Information about updates can be obtained 

• on the World Wide Web: 

http : //crunch . ikp . physik . th-darmstadt . de/nlc/beam . html 

• by internet FTP: 

host: crunch . ikp .physik. th-darmstadt . de 

user: anonymous 

password: your email address 

directory: pub/ohl/circe 

• from mailing lists: 

circe-announceOcrunch . ikp . physik . th-darmstadt . de 
circe-bugsOcrunch . ikp . physik . th-darmstadt . de 
circe-discussOcrunch. ikp. physik. th-darmstadt .de 

Subscriptions are available from 

maj ordomoOcrunch . ikp . physik . th-darmstadt . de 

Contributions of results from other simulation programs and updated accelera- 
tor designs are welcome at 

Thorsten . OhlOPhysik . TH-Darmstadt . de 
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